Polymer collapse in the presence of hydrodynamic interactions 
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We investigate numerically the dynamical behaviour of a polymer chain collapsing in a dilute 
solution. The rate of collapse is measured with and without the presence of hydrodynamic inter- 
actions. We find that hydrodynamic interactions both accelerate polymer collapse and alter the 
folding pathway. 
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I. INTRODUCTION 



When a polymer is placed in solution hydrophobic in- 
teractions between monomers and solvent molecules can 
cause it to undergo a collapse transition to a compact 
state [0,D . The statistical physics of the polymer tran- 
sition from the extended to the collapsed state is well 
understood. However the dynamics of the transition and 
in particular the effect on that dynamics of the hydrody- 
namic properties of the solvent remains unclear. There- 
fore in this Letter we investigate numerically the dynam- 
ical behaviour of a polymer chain collapsing in a dilute 
solution. The collapse is measured with and without the 
presence of hydrodynamic interactions thus allowing a 
direct investigation of their effect. We find that hydro- 
dynamics accelerates the polymer collapse. It also alters 
the folding pathway, allowing the folding to occur more 
homogeneously along the polymer chain rather than ini- 
tially at the chain ends. 

With the introduction of Zimm's model |^ it became 
apparent that hydrodynamics play a central role in the 
dynamics of polymers in dilute solution. However, un- 
derstanding such interactions is difficult, analytically be- 
cause they present a complicated many-body problem, 
and numerically because they develop on time-scales long 
compared to the thermal fluctuations of the monomers. 

Theoretical work on the dynamics of polymer collapse 
can be divided into two approaches. Phenomenological 
models balance the driving and dissipative forces to give 
scaling laws They involve assumptions about how 

the collapsed state develops on which there is no con- 
sensus. Several authors have considered models which 
are based on a solution of the Langevin equation . 
Of particular interest is work by Pitard jl^] and by 
Kuznetsov et al who find the inclusion of hydrody- 
namics, modelled by a preaveraged Oseen tensor, speeds 
up the collapse. 

Simulations on polymer collapse P,|T3|-p^ , using 
Monte Carlo or Langevin approaches, have not included 
the hydrodynamic effects of the solvent. Very recent work 
has shown that it is now possible to use molecular dy- 
namics simulations with an explicit solvent to model the 
collapse transition if powerful computational resources 
are available Q. In an interesting recent paper Chang 



and Yethiraj |T^] compared molecular dynamics simu- 
lations of a polymer in a solvent to Brownian dynam- 
ics simulations. They attempted to match the param- 
eters in the two simulations and hence to compare col- 
lapse with and without hydrodynamics. Here we use a 
hybrid approach [0 where the solvent is modelled by 
a Malevanets-Kapral method and the polymer by 
molecular dynamics in our investigation of the hydrody- 
namics of polymer collapse. 



II. SIMULATION DETAILS 

Modelling a dilute polymer solution is a difficult task 
because of the existence of widely differing time scales. 
The dynamical properties of polymers can be dominated 
by hydrodynamic interactions between different parts of 
the polymer chains In contrast with the time scale 
of thermal fluctuations of individual monomers, these in- 
teractions are long-ranged and evolve slowly. Therefore 
it is computationally too expensive to reach hydrody- 
namic time scales using molecular dynamics simulations 
for both the polymer and the solvent molecules. 

To overcome this problem we use a hybrid simulation 
approach where the equations of motion of the polymer 
alone are solved using a molecular dynamics algorithm. 
The solvent is modelled using a mesoscale approach, de- 
veloped by Malevanets and Kapral [Q. This ignores 
the molecular detail of the solvent but preserves its abil- 
ity to transmit hydrodynamic forces. The polymer can 
be thought of as moving within a "hydrodynamic heat 
bath" . 

The polymer chain is modelled by beads connected via 
non-harmonic springs [l9|] with adjacent beads along the 
chain backbone representing an effective Kuhn length of 
the polymer chain. These flnitely extensible springs are 
represented by the FENE potential 
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A Lennard-Jones potential Ig^l which acts between all 
the polymer beads is used to model the excluded volume 
of the monomers and a long range attraction which drives 
polymer collapse 
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We take e = 1.0, a = 1.0, k — 30 and Rq — 2 where 
parameters and results are quoted in reduced Lennard- 
Jones units. 

Newton's equations of motion for the polymer are inte- 
grated using the time reversible velocity Verlet algorithm 
p0| . The molecular dynamics time step is chosen to be 
6t — 0.002 where ts is the interval between solvent col- 
lision steps, defined below. 

The solvent is modelled by a large number N — 131072 
of point-like particles which move in continuous space 
with continuous velocities but discretely in time |l|] . The 
algorithm is separated into two stages. In the first of 
these, a free streaming step, the positions of the solvent 
particles at time t, Xi{t), are updated simultaneously ac- 
cording to 



{t+ts) = Xi (t) + {t)ts 



(3) 



where Vi{t) is the velocity of a particle. 

The second component of the algorithm is a collision 
step which is executed on both solvent particles and poly- 
mer beads. The system is coarse-grained into Li^ unit 
cells of a regular cubic lattice. In this simulation L = 32 
is used. There is no restriction on the total number of sol- 
vent or polymer particles in each cell, although the total 
number of particles is conserved. Multiparticle collisions 
are performed within each individual cell of the coarse- 
grained system by rotating the velocity of each particle 
relative to the centre of mass velocity Vcm{t) of all the 
particles within that cell 



Vi {t + ts 



(t) + R{Vi{t) ~ Vcm{t) ) . 
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i? is a rotation matrix which rotates velocities by 9 
around an axis generated randomly for each cell and at 
each time step. In the present calculations we take 6 = ^. 

Note that the collision step preserves the position of 
the solvent and polymer beads. It transfers momentum 
between the particles within a given cell while conserv- 
ing the total momentum of these particles. Because both 
momentum and energy are conserved locally the thermo- 
hydrodynamic equations of motion are captured in the 
continuum limit |p^ . Hence hydrodynamic interactions 
can be propagated by the solvent and, because the poly- 
mer beads are involved in the collisions, to the polymer. 
Note, however, that molecular details of the solvent are 
excluded: this allows the hydrodynamic interactions to 
be modelled with minimal computational expense. 

The volume in phase space is invariant under both the 
free streaming and collision steps. Hence the system is 
described by a microcanonical distribution at equilibrium 
I p^ . The initial solvent distribution was generated by 
assigning positions randomly within the system with an 
average density p = 4 particles per unit cell. The ve- 
locities were assigned from a uniform distribution which 



relaxed rapidly {tg < 100) to the equilibrium Maxwell- 
Boltzmann form. 

A particularly useful feature of the Malevanets-Kapral 
algorithm is the ease with which hydrodynamic interac- 
tions can be "turned-off" thus replacing the hydrody- 
namic heat bath by a Brownian (random) heat bath. 
This is achieved by randomly interchanging the veloci- 
ties of all the solvent particles after each collision step, 
thus relaxing the constraint of local momentum conser- 
vation to a global one. Accordingly the velocity corre- 
lations which result in hydrodynamic interactions disap- 
pear from the fluid. Running simulations with the same 
initial conditions and parameter values, but with hydro- 
dynamics present or absent, greatly faciliates pinpointing 
the effect of the hydrodynamic interactions. 



III. RESULTS 

Our aim is to study the dynamics of polymer collapse 
from an extended to a compact state [0. The collapse 
trasnsition is driven by the attractive Lennard-Jones in- 
teractions between the polymer beads. The transition is 
rounded and shifted by the finite chain length, and takes 
place at fc^T 1.8 for chains of length Np = 100. 

Initial polymer configurations were chosen at random 
from long runs on a polymer chain at equilibrium in a 
solvent with ksT = 4. Each of these extended con- 
figurations were placed in a solvent at equilibrium at 
ksT = 0.8. This value ensures collapse but prevents the 
chain collapsing so quickly that hydrodynamic interac- 
tions do not have sufficient time to develop. The rate of 
polymer collapse was measured by monitoring the varia- 
tion of the radius of gyration Rg of the chain with time 
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A typical numerical result is shown in Figure 1 for a chain 
of length Np = 100 with and without hydrodynamics. 

The collapse time was estimated as the time when the 
equilibrium radius of gyration was first attained. Results 
for T are shown for chains of varying lengths in Table 1. 
T was averaged over 20 initial configurations for Np — 40, 
60 and 100, and 5 for Np — 200, and the collapse time 
for hydrodynamic (th) and Brownian (tb) heat baths 
compared. Variations in collapse times between differ- 
ent runs are large as expected. However, it is strikingly 
apparent that hydrodynamic interactions speed up the 
rate of collapse by a factor ^ 2 for each polymer length. 
For all but 6 of the 65 runs and for all the Np — 100 
and 200 runs the collapse was faster with hydrodynamics 
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switched on. The results in Table 1 also illustrate that 
the collapse time increases with chain length as expected 



1.6±0.2 
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). Values predicted in the 



literature for the exponents relating the collapse time to 
chain length vary widely. It is far from obvious that 
these exponents are universal: they appear to depend on 
quench depth and the details of the model More- 
over the final radius of gyration was recorded for both 
the Brownian and hydrodynamic collapses and found to 
agree as the equilibrium polymer properties should be 
independent of the nature of the heat bath. 

Figure 2 compares typical collapse pathways with and 
without hydrodynamics for Np = 100 and a final temper- 
ature ksT = 0.8. In this and the majority of other runs a 
qualitative difference was observed in the collapse mech- 
anism. In the Brownian solvent the ends of the polymer 
tend to collapse first forming a dumbbell shape. With hy- 
drodynamics the collapse takes place more evenly along 
the chain. We speculate that this occurs due to 'slip- 
streaming'. The movement of beads towards each other, 
which is initiated by the attractive Lennard- Jones poten- 
tial, is enhanced by hydrodynamic interactions. Once a 
bead starts moving in a given direction and locally drags 
fluid with it, it is easier for neighbouring beads to move 
in the same direction. Similar behaviour albeit on much 
larger length scales has been observed experimentally by 
Bartlett et al |^l|] for a pair of colloidal particles in solu- 
tion. 



IV. DISCUSSION 

We have shown that it is possible to measure the dy- 
namics of the collapse transition of a polymer in a solvent 
using a hybrid mesoscale/molecular dynamics algorithm. 
Hence it has been possible to show that hydrodynamic 
interactions speed up the collapse of the polymer chain. 
The hydrodynamics alters the collapse pathway, allow- 
ing the folding to occur more homogeneously along the 
chain, rather than initially at the chain ends. 

Qualitative features of the collapse pathways are in 
agreement with those reported by Chang and Yethiraj 
who recently compared molecular dynamics simu- 
lations of a polymer in a solvent, which include hydro- 
dynamics, to Brownian dynamics simulations, which did 
not. These authors also found that the Brownian simu- 
lations could become trapped in a metastable free energy 
minimum. It would be interesting to address whether this 
is due to deeper quenches than those considered here or 
to the different simulation approaches employed. 

The solvent is modelled using the Malevanets-Kapral 
algorithm which sustains hydrodynamic modes but does 
not include molecular interactions. Hence we caution 
that we cannot investigate the late stages of collapse 
where trapping of the water molecules in the polymer 



chain may be important. Moreover collapse is driven by 
attractive Lennard-Jones interactions between the poly- 
mer beads rather than repulsive monomer-solvent inter- 
actions. It would be of interest to include the solvent 
near the polymer in the more realistic molecular dynam- 
ics updating and the possibility of doing this is under 
investigation. 

Finally it may be of interest to consider whether the 
early stages of the folding of proteins in aqueous solution 
might be influenced by the hydrodynamic interactions of 
the solvent increasing the cooperativity of the collapse. 
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FIG. 1. Variation of the radius of gyration with time for a collapsing polymer chain in a hydrodynamic ( — ) or Brownian 
(• • •) heat bath. 
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FIG. 2. A comparison of the pathways for polymer collapse with and without hydrodynamics. 





{th} 


(rB) 






40 


290 ± 92 


401 ±118 


2.02 ± 0.05 


2.02 ± 0.04 


60 


490 ± 79 


671 ± 170 


2.23 ± 0.02 


2.23 ± 0.01 


100 


858 ± 113 


1571 ± 252 


2.58 ± 0.01 


2.58 ± 0.02 


200 


2050 ± 25:! 


5500 = :!<)() 


:!.15 ± 0.00 


± 0.02 



TABLE I. Averaged collapse time in units of the solvent time step of a polymer chain of length A'p with {{th)) and without 
{{tb}) hydrodynamics for a final temperature fcsT = 0.8. The radius of gyration after collapse in reduced Lennard-Jones units 
with ({Rgu)) and without {{Rgg}) hydrodynamics is also listed. 
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